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Abstract 

The recombination of hydrogen in the interstellar medium, taking place on surfaces of micro- 
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£vq . scopic dust grains, is an essential process in the evolution of chemical complexity in interstellar 



clouds. Molecular hydrogen plays an important role in absorbing the heat that emerges during 

o ' 

gravitational collapse, thus enabling the formation of structure in the universe. The H2 formation 



process has been studied theoretically, and in recent years also by laboratory experiments. The 
experimental results were analyzed using a rate equation model. The parameters of the surface, 
that are relevant to H2 formation, were obtained and used in order to calculate the recombination 
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rate under interstellar conditions. However, it turned out that due to the microscopic size of the 
dust grains and the low density of H atoms, the rate equations may not always apply. A master 

*t2 ! 

equation approach that provides a good description of the H2 formation process was proposed. It 
takes into account both the discrete nature of the H atoms and the fluctuations in the number 
of atoms on a grain. In this paper we present a comprehensive analysis of the H2 formation pro- 



cess, under steady state conditions, using an exact solution of the master equation. This solution 
provides an exact result for the hydrogen recombination rate and its dependence on the flux, the 



surface temperature and the grain size. The results are compared with those obtained from the 
rate equations. The relevant length scales in the problem are identified and the parameter space 
is divided into two domains. One domain, characterized by first order kinetics, exhibits high effi- 
ciency of H2 formation. In the other domain, characterized by second order kinetics, the efficiency 
of H2 formation is low. In each of these domains we identify the range of parameters in which, due 
to the small size of the grains, the rate equations do not account correctly for the recombination 
rate and the master equation is needed. 

PACS numbers: PACS:05.10.-a,82.65.+r,98.58.-w 
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I. INTRODUCTION 



The recombination of hydrogen on the surfaces of microscopic dust grains in the inter- 
stellar medium has attracted much interest in recent years. This process is essential since 
gas-phase reactions cannot account for the abundance of H 2 in interstellar clouds ^ |3|, [4] . 
Theoretical |5], ||, |7], |8], [IT, |12[ and experimental [|13|, [14], [15|, [16], [17], |18[| techniques 

have been used in order to evaluate the rate of H 2 formation on relevant dust materials 
under interstellar conditions. Quantum mechanical calculations and molecular dynamics 
simulations were performed, in attempt to identify the diffusion and reaction rates on the 
surfaces of various astrophysically relevant materials 0, [| [|, [], [7], [T2], |2"0fl . Experimental 
results were also obtained for the activation energies of the relevant diffusion and desorption 
processes on various surfaces fl3], [14], |15 |. 

Rate equations are an essential tool in the modeling of chemical reactions in the interstel- 
lar medium |2l], p2l 



H, H 1111,13,11 HI S IS HI- Chemical models based on the rate 
equation approach take into account a large number of reactions in the gas phase as well as 
reactions that take place on the surfaces of dust grains. Rate equations were recently used in 
order to analyze the results of laboratory experiments on H 2 formation on dust-analog sur- 
faces [33|. The analysis provided the surface parameters that are essential for the evaluation 
of the H 2 formation rate on dust grains in the interstellar medium. The rate equations used 
in Ref. |33| describe the diffusion, reaction and desorption processes on the surface. They 
provide the time evolution of the average densities of atoms and molecules on the surface, 
while fluctuations are neglected. Such rate equations are expected to provide good results 
for macroscopic surfaces. However, it turns out that they may be unsuitable for the study 
of hydrogen recombination in the interstellar medium due to the small grain size and low 



flux [34], |35|, [36], |37], |38[ . Under these conditions the number of H atoms on the surface of a 
grain may be very small and fluctuations are expected to be significant. Attempts to resolve 
this difficulty included the use of modified rate equations in which parameters are changed 



to account for the finite grain size J34] , p5| , pq , p7| , pqj . Monte Carlo methods were also used, 
to simulate the surface diffusion and reaction processes on small grains [^, £55). The class 
of Monte Carlo methods that are suitable for such activated processes on the surface, away 
from thermal equilibrium, are the continuous time or kinetic Monte Carlo techniques 



In these simulations, at each time step, the next move is picked with a probability propor- 



2 



tional to its rate. The elapsed time is given, according to the theory of stochastic processes, 
by the inverse of the sum of the rates of all processes that could have occurred at that 
time. The kinetic Monte Carlo approach can be directly related to the underlying master 
equation that describes the time evolution of the probabilities of all the microscopic states 
of the system [I2"|. Monte Carlo simulations typically require large computational resources. 
For example, the calculation of averages for quantities such as the H2 formation rate is done 
by collecting large amounts of statistical information |39|, f40[. To study the chemistry of 
interstellar clouds one needs a model that couples the gas-phase and grain-surface reactions. 
It was found that the use of Monte Carlo methods in this context is impractical, while rate 
equations for surface chemistry cannot account correctly for reaction rates on small grains. 
Recently it was shown that the H 2 formation process on small grains can be described 



by a master equation approach |43| , |44}| . The master equation takes into account both the 
discrete nature of the H atoms as well as the fluctuations. Its dynamical variables are the 
probabilities P(n) that there are n atoms on a grain at time t. The time derivatives P{n), 
n = 0,1,2,... are expressed in terms of the adsorption, reaction and desorption terms. 
The master equation provides the time evolution of the probabilities P{n), from which the 
recombination rate can be calculated. It was used in conjunction with surface parameters 
obtained experimentally and computationally, to explore the hydrogen recombination and 
other chemical reactions on small grains under interstellar conditions [^, . 

In this paper we present a comprehensive analysis of the H2 formation process under 
steady state conditions as a function of the physical parameters. These parameters include 
the flux of H atoms, the grain size and surface temperature. They also include properties of 
the surface, namely the density of adsorption sites as well as the activation energies for H 
diffusion and desorption. The analysis is based on an exact analytical solution of the master 
equation. In this solution, the steady state distribution P(n) is expressed in terms of two 
dimensionless quantities, composed of the physical parameters mentioned above. Using this 
solution we identify the relevant length scales in the problem. The parameter space is then 
divided into two domains. One domain, characterized by first order kinetics (namely, the 
H 2 formation rate is linearly proportional to the flux of H atoms), exhibits high efficiency 
of H 2 formation. In the other domain, characterized by second order kinetics (namely, the 
H 2 formation rate is proportional to the square of the incoming flux), the efficiency of H 2 
formation is low. In each of these domains we identify the range of parameters in which, 
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due to the small size of the grains, the rate equations do not apply and the master equation 
is needed. 

The paper is organized as follows. The rate equation model is described in Sec. IIA and 
analyzed in Sec. IIB. The master equation is introduced in Sec. IIIA. The exact solution of 
the master equation is given in Sec. IIIB. The analysis of the H 2 formation process on small 
grains, that is based on this solution, is presented in Sec. IIIC. The results are discussed in 
Sec. IV and summarized in Sec. V. 



II. H 2 FORMATION ON MACROSCOPIC SURFACES 



A. The Rate Equation Model 

Consider a flux of H atoms that are irradiated and stick to a macroscopic surface. The 
atoms perform hops as random walkers on the surface and recombine into H2 molecules when 
they encounter one another. Let p(t) (in monolayers [ML]) be the coverage of H atoms on 
the surface at time t. Its time dependence is described by the following rate equation: 

f t =f-Wp-2ap\ (1) 

The first term on the right hand side of Eq. ([!]) represents the flux of H atoms. The parameter 
/ represents the effective flux of atoms (in units of ML s~ l ), namely, the (temperature 
dependent) sticking coefficient £(T) is absorbed into / [the Langmuir-Hinshelwood (LH) 
rejection process [IB| of atoms deposited on top of already adsorbed atoms is neglected here 



since it is assumed that the coverage is low]. The second term in Eq. ([]]) represents the 
desorption of H atoms from the surface. The desorption coefficient is 



W = v ■ exp(-E 1 /k B T) (2) 

where v is the attempt rate (standardly taken to be 10 12 s _1 ), E\ is the activation energy 
barrier for desorption of an H atom and T is the temperature. The third term in Eq. (|IJ) 
accounts for the depletion of the H population on the surface due to recombination into H 2 
molecules, where 



a = v • exp(— Eo/ksT) 



(3) 



is the hopping rate of H atoms on the surface and Eq is the activation energy barrier for 
hopping. Here we assume that diffusion occurs only by thermal hopping, in agreement with 
experimental results [Q. We also assume that there is no energy barrier for recombination. 
The H 2 production rate r (ML s _1 ) is given by 



r = a p\ (4) 

For simplicity we assume here that all the H 2 molecules are desorbed from the surface upon 
formation. Even if on real surfaces some of the molecules may remain on the surface and 
desorb thermally later, under steady state conditions, at low coverage, this will not affect 
the recombination rate. 

A (more complete) model, based on Eq. (|1|) was used |53] to analyze the results of 



temperature programmed desorption (TPD) experiments (|TJ, [TJ], |T5|] . The values of Eq 
and Ei (as well as two additional parameters) that best fit the experimental results were 
obtained. The steady state behavior under astrophysically relevant conditions was then 
studied and the recombination efficiency 

" = 775 (5) 

was calculated in an astrophysically relevant range of flux and temperature. Note that 
under steady state conditions rj is limited to the range < rj < 1. It was found that the 
recombination efficiency is highly temperature dependent. It exhibits a narrow window of 
high efficiency along the temperature axis, which slowly shifts to higher temperatures as the 



flux is increased 33 . 



B. Analysis and Results 

Consider a macroscopic surface exposed to a constant flux of H atoms, as described by 
Eq. ([[]). Under steady state conditions dp/dt = and the coverage is given by 



1 (W 
4 \ a 



Xf/W) 

V 1 



1 + 



(6) 



The recombination efficiency rj = 2(a/f)p 2 takes the form 
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''-8(7Mf 1 + ^ 1 + 8 (HW • (7) 

To obtain a better understanding of the recombination process we will try to identify the 
length and time scales involved. Consider the ratio a/W between the hopping rate and the 
desorption rate. This ratio is (on average) the number of hops an H atom visits before it 
desorbs (neglecting recombination). It is also (up to a constant multiplicative factor of order 
unity and a logarithmic correction) the number of sites that the atom visits before it desorbs 
pE7| . We will denote this number by 



Svisit = a/W. (8) 

Consider the ratio W/f between the desorption rate and the flux. Neglecting the recom- 
bination term we obtain W / f = 1/p, namely this is approximately the average number of 
vacant sites around each H atom. We denote it by 



Vacant = W/f. (9) 

The properties of the recombination process can thus be plotted on a two dimensional 
parameter space, namely f /W vs. W/a [Fig. |T[. The coverage p depends on these two 
parameters, while 77 depends only on the ratio 

The coverage p and the recombination efficiency 77 as a function of 9 along the curve 

1 c^y, (id 



W \a 

where c = 1/2, are shown in Figs. |2j(a) and 0(b), respectively. The solid lines in Figs. 0(a) 
and 0(b) show the results without LH rejection, while the dashed lines show the results for 
similar conditions but including LH rejection. In these Figures we identify a domain of high 
coverage and high efficiency for f /W > W/a and a domain of low coverage and low efficiency 
for f /W < W/a, separated by the diagonal line in Fig. |I[ We will now analyze the limits, 
deep in each of these domains, in terms of the length scales associated with s visit and s vacan t 
(the length-scales are obtained by dividing s visit and s vacan t by the density of adsorption sites 
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and taking the square root). In the case / /W 3> W/a, namely s v i S i t ^> s vacant , the typical 
number of sites that an atom visits is much larger than the number of vacant sites around 
it. Therefore, it is most likely to find a second H atom and recombine. In the opposite 
limit f /W <C W/a, namely s visit <C s vacant , most atoms visit only vacant sites around their 
initial adsorption site and desorb before having a chance to form molecules. As a result, the 
recombination efficiency rj is low. Evaluating the coverage in both limits we obtain 

w 

^ J V2 V a ■ W " a / 12 N 

X ■ X <^ K 

{ W ' W ^ a ■ 

The H 2 production rate r = ap2 is given by 

X NN K 

(13) 

In the first limit, r is linear in the flux /, thus we denote this limit as first order kinetics. 
In the second case r is proportional to f2 and we thus denote it as second order kinetics. 
The regime of first order kinetics is characterized by high recombination efficiency rj = 1, 
since the desorption of H atoms is negligible. In the second order kinetics, desorption of H 
atoms (before they form molecules) is the dominant process. Therefore, the efficiency rj is 
low. More precisely, the recombination efficiency is 

1 ■ J- > ^ 
V ={ w a (14) 

1 2s. f ■ X <^ YL 

W 2 J ■ W ^ a ■ 

For a given surface it is convenient to plot the coverage and recombination rate as a 




function of the flux and the surface temperature p3| . The boundary between the first and 
second order domains is given by 

*w> - drab? (15) 

with the first order domain for temperatures below this curve and the second order domain 
above it. For parameters on this boundary line the recombination efficiency is r] = 1/2. 

Note that in Eq. (0) the coverage p is not limited to the range < p < 1. In the LH 
kinetics, the coverage is bounded from above due to the rejection of atoms deposited on top 
of other H atoms that are already adsorbed. To include this effect, the flux term / in Eq. ([!]) 
is replaced by f{l — p). Solving the new equation at steady state we obtain 



_ \ {W + f)/a 
77 8 f/(W + f) 



•1 + 



1 1 ? f/(W + f) 
\ + (W + f)/a 



(16) 



which is similar to Eq. (|?p except that now W is replaced by (W + /). The LH rejection 
introduces further constraints on the domain of high efficiency of H2 formation. This is 
due to the fact that at low temperature the coverage approaches unity and newly deposited 
atoms are rejected. For simplicity we will define the high efficiency domain as the set of 
points in parameter space for which rj > 1/2. The boundary of this domain, namely the 
curve in the (W/a, f /W) plane on which 77 = 1/2, is given by 

^-(TTW (17) 

In Fig. H] we show the high efficiency domain (gray area) and the low efficiency domain in 
parameter space for the LH kinetics. For f /W <C 1 the boundary between them coincides 
with the diagonal line that separates the first and second order domain without the LH 
rejection. For any given value of W/a in the range < W/a < 1/4, the high efficiency 
domain is bounded from above and below by f\ ow (W/a)/W < f/W < f np (W/a)/W where 
/up (W/a) and f\ w(W/a) are determined by Eq. fllTD - The two boundaries are related to 
each other according to 

MW/a) W 

W fUW/aY 1 ) 

For a given surface, for which Eq, E\ and v are known, we can draw a diagram for the 
recombination process in the (T, log 10 /) plane (Fig. ^). For a given flux, the temperature 



in the domain of high efficiency (gray area) is now bounded from above by Eq. fll5D and 
from below by 

7low(/) = , E ° ,y (19) 

k B Qnv- m/) 

Therefore, in order for a given surface to exhibit a domain of high efficiency, the condition 
T U p(f) > ^iow(/) must be satisfied. This condition is satisfied if E\ > Eq, namely the 
activation energy for desorption is higher than for diffusion. 

A further constraint on the high efficiency domain may appear due to the existence of H2 
molecules on the surface. These molecules may reject some of the deposited atoms through 



S 



the LH mechanism and thus reduce the effective flux. The H2 molecules on the surface may 

or 



be either molecules that formed on the surface and did not desorb upon formation |33| 
ones that were adsorbed from the gas phase. In both cases the condition for the existence 
of a high efficiency domain is that E 2 < 1E\ — E , where E 2 is the activation energy for 
desorption of H 2 molecules from the surface. Thus, surfaces that do not satisfy this condition 
(particularly surfaces on which H2 molecules are adsorbed more strongly thatn H atoms) are 
not expected to efficiently catalyze the H 2 formation process, unless molecules are desorbed 
upon formation and their density in the gas phase is too low to saturate the adsorption sites 
on the surface. 

For small grains and low flux one may reach the situation in which the average number 
of H atoms on a grain is of order unity or even less. Under these conditions the rate 
equation model, that takes into account only average densities, ignoring the fact that H 
atoms are discrete entities, does not account correctly for the recombination rate. This is 
due to the fact that the recombination process requires at least two atoms on the surface 
and the fluctuations in the number of H atoms on different grains become dominant. A more 
complete description of the recombination process is needed. Such description is provided 
by the master equation presented below. 



III. H 2 FORMATION ON SMALL DUST GRAINS 



A. The Master Equation Model 

We will now consider the formation of H 2 molecules on small dust grains. In this case it is 
more convenient to rescale the parameters such that instead of using quantities per unit area 
- the total amount per grain will be used. The number of H atoms on the grain is denoted 
by n. Its expectation value is given by (n) = S ■ p where S is the number of adsorption sites 
on the grain. The incoming flux of H atoms onto the grain surface is given by F = S ■ f 
(atoms s" 1 ). The desorption rate W remains unchanged. The hopping rate a (hops s _1 ) 
is replaced by A = a/S which is approximately the inverse of the time t s required for an 
atom to visit nearly all the adsorption sites on the grain surface. This is due to the fact 
that in two dimensions the number of distinct sites visited by a random walker is linearly 
proportional to the number of steps, up to a logarithmic correction [[|7j . The H 2 production 
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rate of a single grain is given by R = S ■ r (molecules s 1 ). The rate equation will now take 
the form 

= F- W(n) - 2A(n) 2 . (20) 

Under given flux and surface temperature, for grains that are large enough to hold many 
H atoms, Eq. (p0|) provides a good description of the recombination process. However, for 
small enough grains (n) becomes or order unity and Eq. (|20|) becomes unsuitable, because 
it neglects the fluctuations in the number of atoms on a grain. 

We will now introduce the master equation, that provides a correct description of the 
recombination process even in the limit of small grain sizes and low flux. Consider a grain 
that is exposed to a flux F of H atoms. At any given time the number of H atoms adsorbed 
on the grain may be n = 0, 1, 2, . . . , S. The probability that there are n hydrogen atoms on 
the grain is given by P(n), where 

EP(n) = l. (21) 

n=0 

The master equation provides the time derivatives of these probabilities, P(n), namely the 
gain or loss of the probabilities of the different states. These derivatives are linear in the 
probabilities themselves. The equations include three terms. The first term describes the 
effect of the incoming flux F. The probability P(n) increases when an H atom is adsorbed on 
a grain that already has n — 1 adsorbed atoms [at a rate of FP(n — 1)], and decreases when 
a new atom is adsorbed on a grain with n atoms on it [at a rate of FP(n)}. The adsorption 
process is considered as completely random (Poisson process), and is fully characterized 
by F. The second term includes the effect of desorption. An H atom desorbed from a 
grain with n adsorbed atoms, decreases the probability P(n) [at a rate of nWP(n), where 
the factor n is due to the fact that each of the n atoms can desorb], and increases the 
probability P(n — 1) at the same rate. The third term describes the effect of recombination. 
The production of one molecule reduces the number of adsorbed atoms from n to n — 2. 
For a given pair of H atoms, the recombination rate is proportional to the sweeping rate A 
multiplied by 2 since both atoms are mobile simultaneously. This rate is multiplied by the 
number of possible pairs of atoms, namely n{n — l)/2. The master equation exhibits the 



Markov property, namely, no memory effects are included [42 . This property emerges from 
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the fact that the incoming flux keeps washing out any spatial correlations that may develop 
due to recombination events between adjacent atoms. The master equation thus takes the 
form: 



p(0) = -FP(0) + WP(1) + 2 ■ 1 ■ AP(2) 
P(l) = F [P(0) - P(l)] + W [2P(2) - P(l)] + 3 • 2 • AP(3) 
P(2) = F [P(l) - P(2)] + W [3P(3) - 2P(2)] 
+ A [4 • 3 • P(4) - 2 • 1 • P(2)] 

P(n) = F [P(n - 1) - Pin)} + W [(n + l)P{n + 1) - nP{n)} 
+ A[(n + 2){n + l)P(n + 2) - n(n - l)P(n)] 

P(S) = F [P(S - 1) - P(S)] - SWP(S) - SiS - l)AP(S). (22) 

Note that the equations for P(0) and P(l) do not include all the terms, because at least one 
H atom is required for desorption to occur and at least two for recombination. Similarly, 
the equation for P(S) does not include all the terms since there is no room for more than S 
atoms on the grain surface. The expectation value for the number of H atoms on the grain 
is 

s 

( n ) = ^nP(n). (23) 

n=0 

The rate of formation of H 2 molecules, R (molecules s -1 ), is thus given by 

s 

R = AY,n(n-l)P(n). (24) 

n=2 

The hydrogen recombination efficiency on the grains is given by 

"=<m- (25) 



n 



B. Exact Solution of the Master Equation 



When a grain is maintained at a constant temperature (namely W and A are fixed), and 
is exposed to a constant flux F, the recombination process on its surface approaches a steady 
state. Under steady state conditions the time derivatives on the left hand side of Eq. ( p^) is 
zero. We thus obtain a homogeneous set of coupled linear equations in the variables P{n), 
n — 0, 1, 2, . . . , S. This set can be expressed in a matrix form as 

MP = 0, (26) 

where 



M 



-F 


w 


2A 








F 


—F — W 


2W 


3 • 2A 








F -F 


- 2(W + A) 


3W 


S-(S-1)A 








F 


-F -3(W + 2A) .. 


SW 











F 


. -F-S[W + (S- 



V 

(27) 

and P consists of the components P(n), n = 0, 1, . . . , S. The matrix elements are denoted 
by M(n, m), n,m = 0,1,2, S . The only non-vanishing matrix elements are 



M(n,n) = -[F + nW + n(n- 1)A) 
M(n + l,n) = F 
M(n,n + 1) = (n + l)W 

M(n,n + 2) = (n + 2)(n + l)A (28) 



For a finite grain the set of equations (]22|) is truncated by n < S. However, under interstellar 
conditions we expect low coverage of H atoms on the grain, namely (n) S. Therefore, 
one can impose a cutoff at some value s < S such that P(n) = for n > s and the 
normalization of probabilities now takes the form J2n=o P( n ) = 1- The terms in the matrix 
that represent flow of probabilities between P(n), n < s and P{n), n > s are removed. 
Three of these terms disappear since they are outside the (s + 1) x (s + 1) size matrix that 

12 



we now consider. The term (-F) in the matrix element M(s, s), that represents the addition 
of an H atom to a grain that already includes s atoms is also removed. The modified term 
will be M(s, s) = -s[W + (s- l)A]. 

We proceed by performing linear operations on the rows of the matrix M. Starting from 
the top, we add each element M(0, m), m = 0, . . . , s of the first row to the corresponding 
element M(l,m) in the second row. We then proceed downwards, adding the elements of 
the nth row to the corresponding elements of the (n + l)th row. Each row of the resulting 
matrix includes one diagonal term and two off-diagonal terms of the form: M'(n,n) = —F, 
M'(n,n + 1) = (n + 1){W + nA) and M'(n, n + 2) = (n + 2)(n + I) A (except for the 
last row in which all the elements are zero and the next to the last row in which the 
second off-diagonal element is removed). In order to remove the second off-diagonal terms 
M'(n, n + 2), n = 0, . . . , s — 2 we now perform a second set of operations, this time starting 
at the bottom rows. We first subtract from each element M'(s — 2,m), m = 0, . . . , s of 
the (s — 2)th row, the corresponding element M'(s — 1, m), of the (s — l)th row, multiplied 
by M'(s — 2, s)/M'(s — 1, s). We then proceed in a similar fashion all the way up. Each 
line in the resulting matrix has only one diagonal and one off-diagonal term. The diagonal 
elements take the form M"(n, n) = —F, n — 0, . . . , s — 1 and M"(s, s) = 0. The off-diagonal 
elements will be 

M"(s-l,s) 
M"(s-2,s- 1) 

M"(s-3,s-2) 
M"(s — n, s — n + 1) 



M"(0, 1) 

By combining the operations in Eq. (^) we express the non-vanishing off-diagonal elements 
of the matrix M" as continued fractions. Their reduced form is 



= M'(s - l,s) 

= M-( S -3, S -2)- ^|i^ .M'( S -2, S -2) 
= M'(s — n, s — n + 1) 



= M'(0, 1) - M ' ( ° ,2 | • M"(l, 1) (29) 
v ' ; M"(l,2) v ' ; v ; 
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M " {n ' n + 1) l^-)^ i — (30, 



n 



+ i)Vaf Vf k a jA { w + n + 1) + _ 



for n = 0, . . . , s — 1. The equation M"P = now takes the form of a set of recursion 
equations: 

Using these equations we can express all the probabilities in terms of P(0) according to 

n - lr M"(i,i) 



p(n) = (-i)-p(o) n 

i=0 



n = l,..., a. (32) 



M"(i,i + 1) 

The probability P(0) is then determined by the normalization condition X)* =0 P(n) = 1. 
Eq. (|32|), complemented by the normalization condition, provides an exact solution of the 
master equation under steady state conditions for any finite cutoff s < S. Since the master 
equation is of use when the coverage is very low, the tail of P(n) already decays for some 
n <C S. Therefore, in evaluating P(0) it is sensible to ignore the cutoff at s < S and write 
the infinite sum instead. In this case, the solution of the master equation can be expressed 
in terms of Bessel functions. The connection to the Bessel functions can be obtained from 
the (infinite) continued fraction expression in Eq. (|3"DD, obtained when s — > oo. Using the 
continued fraction expansion of the ratio J v (z) / J v -i(z) in Ref. |49| (Eq. 9.1.73 in page 363) 
and the relation l v {z) j l v -\{z) = —iJ v {iz)j J v -\{iz) we obtain that 



M"{n,n + 1) _ Iw/A+n-i (2\[F~/A 



(n + l)VFA I w/A+n (2y/F/A). 
We thus obtain an expression for P(n) of the form 



(33) 




P(„)= P( ) W- -L A v L-L . (34) 



The normalization factor P(0) can be expressed in terms of Bessel functions, using Ref. |49 
(Eq. 9.6.51 in page 377, with A = a/2), as 



^-Vff. (35) 
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Therefore, 



fW = ^!f 1 g)'''^W, (36) 



/?■! 




in agreement with Ref. |44[], where the solution was obtained using a generating function. 
To examine the effect of the cutoff at s < S, we compared the distributions P(n) obtained 
from Eqs. (|32|) for different values of s, as well as the distribution expressed in terms of 
the Bessel functions for which s — > oo. Under the conditions of low coverage studied here 
we observe a very fast convergence of Eq. (|3^) to Eq. (|36|) as s increases. We have also 
performed direct numerical integration of the master equation and found that the solution 
given by Eq. (^) is stable and the convergence of the integration process is fast. 

Using similar summations, we can now find exact expressions for the first and second 
moments of the distribution P(n). The average number of H atoms on the grain is given by 



2A hv/A-i(2^2F/A 
The rate of formation of H 2 molecules is given by R = A(n(n — 1)), where 



F Vi (^2F/A 



(n(n - 1)) = — 1 ) V * (38) 
and the recombination efficiency is 

Iw/A-i (2 V /2FM) 

The second moment (n2) can be obtained as the sum of the right hand sides of Eqs. ( pTf ) 
and fl38|). The fluctuations in the number of H atoms on different grains can be quantified 
by the standard deviation of the distribution P(n), given by 



a = ^/(n 2 ) - (n) 2 . (40) 

The limits of first and second order kinetics for large grains can now be reproduced from 
Eq. (|39D . The extreme limit of first order kinetics is characterized by negligible desorption, 
namely W/A 1 and for large enough grains this limit also satisfies F/W ^> 1. In this 
limit the indices of the Bessel functions in Eq. (p9| ) W/A ± 1 — > ±1, respectively, and due 
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to the symmetry I_ n (z) = I n (z) of the Bessel functions |49[ (Eq. 9.6.6 in page 375) 77 — >• 1, 
in agreement with the rate equations. The extreme limit of second order kinetics on large 
grains is characterized by F/W <C W/A. Assuming that the geometric mean of F/W and 



W/A satisfies JF/A « 1 we obtain 

2F 1 

In the case of very large grains, namely W/A > 1 we thus obtain rj = 2AF/W2 = 2af /W2 
which is in agreement with the rate equation results. 

Using Eqs. (p3|) and (|22|) to express the time derivative of (n) we obtain 

*& = f-W(n)-2A(n{n-l)). (42) 

This equation resembles the rate equation (^0| ) except for the recombination term in which 
(n)2 was replaced by (n2) — (n). In the limit of small grains, where (n) is small while the 
fluctuations represented by a become dominant, the rate equation becomes unsuitable and 
over-estimates the rate of H 2 production. 



C. Analysis and Results 

Consider the recombination process on a small grain with S adsorption sites under steady 
state conditions. The boundary between the first order and the second order regimes can 
be expressed in terms of F = f ■ S and A = a/S, taking the form F/W = W/A (neglecting 
the LH rejection). In the first order regime F/W ^> W/A, while in the second order regime 
F/W <C W/A. The domains of first and second order kinetics are shown in Fig. ||, separated 
by the diagonal line F/W = W/A. 

The finite size of the grain introduces a third length scale to the problem. We will now 
examine how small the grain should be in order for the recombination efficiency to deviate 
significantly from the rate equation results. To this end we calculate rj as a function of 
the grain size S, using the exact solution of the master equation, presented above. The 
efficiency rj vs. S in the case of first order kinetics is shown in Fig. |6](a). The recombination 
efficiency obtained from the master equation (solid line) coincides with the rate equation 
result (dashed line) for large grains but declines below some grain size. In the first order 
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case, such deviations typically occur only for extremely small grains of a few thousands 
adsorption sites. The average number of atoms on the grain vs. S is shown in Fig. |](b). 
The S axis in Figs. ^| corresponds to the arrow drawn in the first order domain of Fig. |5|. 
Identifying the corresponding symbols, we observe that the significant decline in rj starts 
when F/W < 1, namely when S < W/f. The fluctuations in the number of H atoms on 
a grain can be quantified by the standard deviation a of the distribution P(n), given by 
Eq. (fief ). The standard deviation a, divided by (n), is shown in Fig. |](c) as a function of the 
grain size S. Clearly, as the grain size decreases, the fluctuations become more pronounced. 

Using the notation introduced above, in the domain of first order kinetics the rate equation 
results start to deviate from the correct value of r\ for grain sizes that satisfy 

S < S vacan t < S visit- (43) 

Under these conditions there is typically no atom or only a single H atom on the grain, 
which is unlikely to find a second atom to recombine with. It thus desorbs before it has a 
chance to recombine. 

The efficiency rj in the case of second order kinetics is shown in Fig. 0(a). Again, the 
master equation efficiency (solid line) coincides with the rate equation result (dashed line) 
for large grains but declines below some grain size, typically of a few tens of thousands of 
adsorption sites. The average number of atoms on the grain vs. S is shown in Fig. 0(b). 
The deviations between the rate equations and the master equation results are accompanied 
by large fluctuations in the number of atoms on a grain, as can be seen in Fig. 0(c). The 
S axis in Figs. [7] corresponds to the arrow drawn in the second order domain of Fig. 
Identifying the corresponding symbols, we observe that the significant decline in rj starts 
when W/A < 1, namely S < a/W . Thus, in the domain of second order kinetics, deviations 
between the rate equations and the master equation occur for a range of grain sizes given 
by 

S <^ s v ; s ;t <i s vacan t. (44) 

In this case, the grain surface area is smaller than the area that an atom can scan before it 
desorbs. As a result the atom tends to perform several sweeps of the grain surface visiting 
again the same vacant sites it has already visited. The probability to find a second atom 
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in these sites is much lower than predicted by the rate equation that does not include such 
return visits. The recombination efficiency is thus sharply reduced as the grain size further 
decreases. 

In conclusion, we observe that in both the first and second order kinetics the recombi- 
nation efficiency on a grain (given by the master equation) starts to deviate from the rate 
equation result when the grain size (represented by S) becomes the smallest length scale in 
the problem. In the first order case this happens when S < s vacant while in the second order 
case it happens when S < s v isit- 

The distribution P(n) at three points along the S axis, in the second order domain, 
are shown in Fig. ^|. For a relatively large grain, P(n) exhibits a well defined and nearly 
symmetric peak. For a very small grain it becomes a a monotonically decreasing function, 
dominated by n = 0, 1 and 2. 



IV. DISCUSSION 



The analysis above can be related to the modified rate equations studied in Refs. [3b 



38| . The modification is needed when the grain size becomes smaller than the two length 
scales involved in the recombination process. In the case of first order kinetics this occurs 
when S < W/ f < a/W, while in second order desorption it occurs when S < a/W < W/f. 

Consider the case of first order kinetics. When the grain becomes smaller than W / f the 
typical number of atoms on the grain is smaller than one (even if we consider the depletion of 
the H atoms on the surface due to desorption alone and neglect recombination). Therefore, 
the rate in which atoms find each other on the surface is not determined anymore by the 
hopping rate a, or the corresponding length scale a/W, but by the grain size 5*. Therefore, 
the recombination term in Eq. fl2"0|) should be modified to reflect the change a/W — > S, or 
equivalently A/W — > 1. This is achieved by replacing A by If. 

In the case of second order kinetics, when the grain size becomes smaller than a/W the 
atom is typically able to perform more than a full sweep of the entire grain before it desorbs. 
However, when an atom visits the same sites for the second time, the probability of finding 
another H atom there is greatly reduced. Therefore, the recombination rate is determined 
by S rather than by a/W, requiring the modification a/W — > S, or the replacement of A 
by W in Eq. (pOf), exactly as in the case of first order kinetics. The modified rate equations 
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thus takes the form 



d <">-- W { n)-2 g (i,%) { nf. (45) 



dt \W F 

where 

I A W ) = [A. #<! or f <1 
"V»" W [W: ^>l and f >1. 
The modified rate equations takes into account correctly the length scales involved in the 
recombination process on small grains. It provides results for (n) that are in significantly 
better agreement with the master equation, compared to the unmodified rate equations. 
It also improves the results for the recombination efficiency However, the modified rate 
equations still involve only the average number of atoms on a grain and do not take into 
account the discreteness of the H atoms and the fluctuations, quantified by the second 
moment in Eq. fl42l) . Since the fluctuations dominate the recombination process on small 
grains, the results of the modified rate equations for rj are not expected to coincide with 
those of the master equation, but only to approximate them better than the ordinary rate 
equations. We observe that the deviations between the rate equation results and the correct 
results obtained from the master equation are significant mostly in the second order domain 
[Fig. [7] . In the domain of first order kinetics such deviations occur only for extremely small 
grains, that may be physically irrelevant [Fig. |6|. 

While the results of the calculations above were presented using dimensionless parameters, 
the actual calculations in Figs. |6| and [7| were done for physically relevant parameters. In both 
of them we used the parameters of the amorphous carbon sample, measured experimentally 



in Refs. [y^, |33[ . On this sample the activation energies for H diffusion and desorption were 
found to be E = 44.0 meV and E\ = 56.7 meV, respectively. The density of adsorption 
sites on the amorphous carbon surface was found to be Sdens = 5 x 10 13 (sites cm -2 ). The 
results for first order kinetics (Fig. |6|) were obtained for T = 17 K and / = 5 x 10~ 8 ML s -1 . 
The results for second order kinetics (Fig. [5]) were obtained for T = 18 K and / = 3.4 x 10 -9 
ML s -1 . The connection with the density and temperature of the hydrogen in the gas phase 
is made through / = p gas v gas /4:Sdens where p gas (atoms cm -3 ) is the density of H atoms in the 
gas phase, v gas is the typical velocity of these atoms and the factor or 4 in the denominator 
is the ratio between the surface area and the cross section for a spherical grain The 
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number of adsorption sites on a spherical grain of diameter d is given by 

S = **(fj "Sdens, (47) 

and the flux F = n(d/2) 2 p gas v gas . Observations indicate that the population of carbonaceous 
and silicate grains in interstellar clouds exhibits a broad distribution of grain sizes, roughly 
in the range 10~ 6 cm < d < 10~ 4 cm [50, |5"T|, BB| 



V. SUMMARY 



We have studied the process of hydrogen recombination on small grains in the interstellar 
medium using the master equation approach. An exact solution of the master equation under 
steady state conditions was presented. This solution provides the probability distribution 
P(n) of having n atoms on the grain as a function of the grain size, flux, temperature and 
the parameters of the surface. From this distribution one can obtain an exact expression for 
the hydrogen recombination rate on the grain surface. The results were compared to those 
obtained from the rate equations that describe the recombination process on macroscopic 
surfaces. In the case of a macroscopically large surface, two length scales are identified. One 
length scale is related to the average number of vacant sites around each adsorbed atom, 
given by s vacant . The other length scale is given by the average number of sites, s visit , that 
an adsorbed atom visits before it desobs (neglecting recombination). The relation between 
these two (length) scales determines the properties of the recombination process, dividing 
the parameter space into two domains. The domain s visit > s vacant is characterized by first 
order kinetics with high recombination efficiency, while the domain s v ; S i t < s vacant exhibits 
second order kinetics with low recombination efficiency. In both domains the finite size of 
the grain enters as a significant factor (requiring the use of the master equation rather than 
the rate equations) when it becomes the smallest (length) scale in the system. In the domain 
of first order kinetics the grain size becomes the smallest (length) scale when S < s vacant 
while in the second order it occurs when S < s^t- 
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FIG. 1: (a) The phase diagram of hydrogen recombination on macroscopic surfaces in the 
(W/a, f/W) plane, as described by Eq. @. The parameter space is divided into two domains: the 
first order domain, above the diagonal and the second order domain below it; 

FIG. 2: (a) The coverage p given by Eq. @ as a function of the angle 6, given by Eq. ( |i"oD (see 
Fig. H), without LH rejection (solid line) and with LH rejection (dashed line). The results shown 
are for points along the curve given by Eq. (|i~T|). (b) The recombination efficiency r\ [Eq. (|7|)] vs. 6, 
without LH rejection (solid line) and with LH rejection (dashed line), along the same curve. Note 
that without LH rejection r\ depends only on 9, namely, it is independent of the specific curve. 

FIG. 4: The phase diagram of the hydrogen recombination process on macroscopic surfaces in 
the (T, In/) plane under the LH kinetics. The strip of high efficiency is now bounded from both 
sides, namely for any given flux there is a range of temperatures that exhibits high recombination 
efficiency. 

FIG. 5: The phase diagram of the hydrogen recombination process on grains without LH rejection, 
in terms of the dimensionless parameters W/A and F/W. The diagonal line separates between 
the domain of first order (on the left) and second order recombination (on the right), respectively. 
The (unit) square near the origin is the domain in which the hydrogen recombination efficiency on 
grains (obtained from the master equation) deviates significantly from the rate equation results. In 
the case of first order processes, the condition for such deviation is F/W < 1, while in second order 
processes the condition is W/A < 1. Also included are two axes, in the first and the second order 
domains, that under given physical conditions (fixed values of / and T) represent the variation in 
the grain size. The behavior of rj, (n) and cr/(n) along these axes is shown in Figs. || and for the 
first and second order cases, respectively. 

FIG. 3: The phase diagram of hydrogen recombination on macroscopic surfaces under the LH 
kinetics. The high efficiency domain is bounded to the gray area, on the left side of the curve given 
by Eq. (0). 



24 



FIG. 6: The hydrogen recombination efficiency r\ (a) and the average number of atoms on the 
grain (n) (b) for first order kinetics as a function of the number of adsorption sites, S, on the grain. 
The solid line (with the symbols on it) shows the results obtained from the master equation, and 
the dashed line shows the rate equation results. The symbols correspond to those in Fig. ||, and 
indicate that the deviations between the master equation and the rate equations become significant 
below F/W = 1. (c) The standard deviation a of the distribution P[n) (normalized by (n)). It 
increases sharply as the grain size is reduced entering the range in which F/W < 1, indicating that 
in this range fluctuations play an important role. 



FIG. 7: The hydrogen recombination efficiency r] (a) and the average number of atoms on the 
grain (n) (b) for second order kinetics as a function of the number of adsorption sites, S, on the 
grain. The solid line (with the symbols on it) show the results obtained from the master equation, 
and the dashed line show the rate equation results. The symbols correspond to those in Fig. ||, and 
indicate that the deviation between the master equation and the rate equations become significant 
below W/A = 1. (c) The standard deviation a of the distribution P[n) (normalized by (n)). It 
increases sharply as the grain size is reduced entering the range in which W/A < 1, indicating that 
in this range fluctuations play an important role. 



FIG. 8: The distribution P(n) at three points along the arrow in the second order domain of Fig. 
||, namely for a small grain with S = 2513 (solid line), and medium size grain with S = 56705 
(dashed line) and for a large grain with S = 141372. 
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